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Abstract 

We present a code that calculates the relic density of the lightest supersymmetric particle (LSP) 
in the minimal supersymmetric standard model. All tree-level processes for the annihilation of 
the LSP are included as well as all possible coannihilation processes with neutralinos, charginos, 
sleptons, squarks and gluinos. In all we have included over 2800 processes not counting charged 
conjugate states. The cross-sections extracted from CompHEP are calculated exactly. Relativis¬ 
tic formulae for the thermal average are used and care is taken to handle poles and thresh¬ 
olds by adopting specific integration routines. The input parameters can be either the soft 
SUSY parameters in a general MSSM or the five parameters of a SUGRA model. A link with 
ISASUGRA/ISAJET allows to calculate the parameters of the MSSM at the weak scale for an input 
at the GUT scale. The Higgs masses are calculated with FeynHiggsFast. Routines calculating 
[g — 2)^ and b ^ are also included. 
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1 Introduction 


One of the very attractive arguments in favour of supersymmetry(SUSY) is that it provides 
a natural solution to the dark matter problem. In R-parity conserving supersymmetric 
models there exists a neutral stable particle, the lightest supersymmetric particle(LSP), 
which could constitute the cold dark matter in the universe. As there are strong con¬ 
straints on stable charged particles, only if the LSP is a neutralino or a sneutrino could 
SUSY provide a suitable dark matter candidate. Cosmological constraints however as 
well as direct searches disfavour the sneutrino |j^ leaving the neutralino as the preferred 
candidate for dark matter. 


The contribution of neutralinos to the relic density is however very model dependent 
and varies by several orders of magnitude over the whole allowed parameter space of 
the MSSM. The relic density then imposes stringent constraints on the parameters of 
the MSSM often favouring solutions with light supersymmetric particles. It has been 
known for some time that coannihilation processes where the LSP interacts with only 
slightly heavier sparticles can significantly reduce the estimate of the relic density, leading 
to acceptable values even with a rather heavy sparticle spectrum. In principle these 
coannihilations can occur with any supersymmetric particleP|. The importance of the 
coannihilation channels were emphasized before both for gauginos P, sleptons^, 
or stops P). Here, ALL channels will be included. In some regions of the parameter 
space the neutralinos annihilate so fast that they cannot constitute the only source of 
dark matter. This happens for example when the masses are such that the neutralinos 
can annihilate through a s-channel Higgs resonance p, or a s-channel Z resonance ||TT 


In this context the inclusion of the relic density constraints, by giving a handle on the 
supersymmetric particle spectrum, has important consequences both for studies of SUSY 
at colliders and in astroparticle experiments. 

There exist many calculations of the relic density in supersymmetry, using various ap¬ 
proximations both in the evaluation of cross-sections and in solving the density equation 

T6[| . Among these, Neutdriver[^]and DarkSusy||13|| are publicly 


12|,|13],[14|,|15 


available. Our purpose is to provide a tool that evaluates with high accuracy the annihi¬ 
lation cross-sections even in regions near poles and thresholds, that is both flexible and 
upgradable and that goes beyond DarkSusy as far as the calculation of the relic density 
is concerned. The main characteristics of this program, called micrOMEGAs, are 


• Complete tree-level matrix elements for all subprocesses 

• Includes all coannihilation channels with neutralinos, charginos, sleptons, squarks 
and gluinos. 



































Loop-corrected Higgs masses and widths 


• Speed of calculation 


All calculations of cross-sections are based on CompHEP[|^, an automatic program for the 
evaluation of tree-level Feynman diagrams. We follow the method proposed by []^| for 
the calculation of the relic density with its generalization to the case of coannihilations 
0. We still rely on approximations for the solution of the relic density equations and the 
determination of the freeze-out temperature, since this allows to significantly increase the 
speed of the program and proves to be very useful when scanning over a large param¬ 
eter space. This, together with the fact that we have included sfermion coannihilation 
channels as well as one-loop corrections to the Higgs width constitute the main differ¬ 
ences with DarkSUSY. Although we will generally assume that the neutralino is the LSP, 
micrOMEGAs can be used to compute the relic density with any supersymmetric particle 
as the LSP, in particular the sneutrino. This is because all (co-)annihilation of any pairs of 
supersymmetric particles into any pairs of standard model or Higgs particles are included. 


The program for the relic density calculation described here is contained in a package 
that lets the user choose between weak scale parameters or parameters of SUGRA models 


as input parameters. The latter is achieved through a link with ISASUGRA/IsajetlllQ 
The calculation of the Higgs masses are done with FeynHiggsFast pO[| . Loop QCD correc¬ 
tions to the Higgs partial widths into fermions are extracted from HDECAY|^I|. In addition 
we provide subroutines that calculate various constraints on the MSSM parameters: di¬ 
rect limits from colliders, Ap, 6 —> sy and {g — 2)^. All these constraints can be updated 
or replaced easily. 


The total number of processes which can contribute to the relic density exceeds 2800. 
However, due to a strong Boltzmann suppression factor, only processes with SUSY par¬ 
ticles close in mass to the LSP are relevant for the calculation. Therefore in most cases 
only a small fraction of the available processes are needed. In principle, compilation of 
the full set of subprocesses is possible, but such a program would be huge and could 
not be distributed easily. To avoid this problem, we include in our package the program 
CompHEPj^ which generates, while running, the subprocesses needed for a given set of 
MSSM parameters. The generated code is linked during the run to the main program 
and executed. The corresponding “shared” library is stored on the user disk space and is 
acessible for all subsequent calls, thus each process is generated and compiled only once. 
Such approach can be realized only on Unix platforms which support dynamic linking. 


The paper is organized as follows. After we summarize the important equations for 
the calculation of the relic density, we give a short description of the parameters of the 
supersymmetric model. A description of the package follows. Finally we present some 






















results and comparisons with another program in the public domain, DarkSusy. 


2 Calculation of the relic density 


The relic density at present in units of the critical density, pcut = SH^/SttG, can be 
expressed as 

m^on^o m^osoYo 


^ Pcrit Pcrit 


( 2 . 1 ) 


where m^o is the mass of the lightest neutralino, H = lOOh kms“^Mpc“^ is the Hubble 
constant and G Newton constant. The entropy conservation assumption makes it easier 
to work with the abundance Y rather than the number density, n^o. sq = s{Tq) dehnes 
today’s entropy to be evaluated at Tq = 2.726K, the temperature of the microwave 
background, with 

27r2 


s(r) = h,n{T)—T 


( 2 . 2 ) 


where hes is a function that depends slowly on the temperature T |^. The code in fact 
returns the relic density 


n^oh^ = 2.755 X 10 ®-^Ho 
GeV 


(2.3) 


One thus needs to hnd Yq = Y(T = Tq). The most complete formulae for the calcula¬ 
tion of Y (T) were presented in [^, ^ and we will follow their approach rather closely. 

The evolution equation of Y is 


dY 

~dT 


45G 


< <7.1, > (V^ - Y^) 


(2.4) 


g^.{T) is a degrees of freedom parameter derived from the thermodynamics describing 
the state of the universe [^, ^ and Yeq = Yeq{T) represents the thermal equilibrium 
abundance 

where we sum over all supersymmetric particles i with mass m* and Qi degrees of freedom. 
Kn is the modihed Bessel function of the second kind of order n . Note that Yeq falls 


rather rapidly as the temperature decreases. < an > is the relativistic thermally averaged 
annihilation cross-section 


/ _ds^Ki{^/T)pl-aij{s) 
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( 2 . 6 ) 


























where cr^ is the total cross section for annihilation of a pair of supersymmetric particles 
into some Standard Model particles, and pij is the momentum of the incoming particles 
in their center-of-mass frame. 

1 

PP = 2 

The summation is over all supersymmetric particles. Integrating Eq. |2.4| from T = oo to 
T = Tq would lead Yq- 


(s - {rrii + mj)‘^){s - (m* - 


(2.7) 


2.1 Preeze-out and approximate solution 


Although one can solve for Y numerically, the procedure is extremely time consuming 
especially when scanning over a large parameter space and when we include a great 
number of processes. It is therefore important to seek as good an approximation as 
possible to speed up the code. We will follow the usual procedure of dehning a freeze-out 


temperature At high T, the LSP are very close to equilibrium and thus Y ~ Wg- 

This will hold until freeze-out, where Y will be almost constant whereas Wg will decrease 
signihcantly. At high T, one can make the approximation that d{Y — Yf.q)/dT is negligible. 
The freeze-out temperature Tj can be dehned from Yf = Y(Tf) = (1 -|- 6)Yeq(Tf) with 6 
some (small) constant number, Tf can then be extracted by solving 


dln{Y, 


eqj 


dT 


45G 


< av > Yeg5{6 + 2) 


( 2 . 8 ) 


In the second regime, where Y S> Y^q, one can neglect Yj‘^ completely. One £nds|]T8 


r(0) 




(2.9) 


In our code, we in fact integrate from T = 0 rather than T = Tq. This above approach 

]. However one usually 


to hnding an approximation was used in numerous papers p5 


omits the 1/Yf term but choose 5 = .5. The authors of ref. state that Yq as given 
in ( p.9|) with 6 = 1.5 agrees better with the precise numerical solution. We checked 
this statement for a wide set of MSSM parameters and reached the same conclusion. 
The approximation agrees with the exact result of a full numerical solution within 2% [|. 
Furthermore we find that the solution (p2.9|) does not depend signihcantly on 6, varying 6 
in the range 1 < 5 < 2 only affects the result by at most 1%. 


It is convenient to express the freeze-out temperature T/ in term of 


^/ = 




T. 


f 


( 2 . 10 ) 


^ The widely used approximation without the 1/1/ term in Eq. 2.9 and S = 0.5 sometimes gives an 
error of about 10%. 



































Typically Xf k, 25, this provides a good initial value for the iterative solution of Eq. |2.8| . 
Only a few iterations are necessary to hnd a solution that satishes Eq. p.8| with 6 = 1.5±.2. 


Typical freeze-out temperatures vary between IGeV and lOGeV, in this temperature 

Thus, if one needs an 


range, the he// and functions can vary by about 20% 


accuracy better than 10%, one cannot use a constant value for these functions as it is 
done in some papers 0. In our program we use the numerical tables of the DarkSUSY 


package [|I^ for a precise evaluation of he// and g^. 


2.2 Numerical integration and summation 


In order to hnd T/ by solving Eq. |2.8| , we have to evaluate several integrals and perform 
a summation over different annihilation channels. In the evaluation of the thermally 
averaged cross-section we have included all two-body subprocesses involving two LSP’s, 
the LSP and a co-annihilating SUSY particle (all the particles of the MSSM) as well as 
all subprocesses involving two coannihilating SUSY particles. The hnal states include 
all possible standard model and Higgs particles that contribute to a given process at 
tree-level. The total number of processes exceeds 2800 not including charged conjugate 
processes. In practice processes involving the heavier SUSY particles contribute only 
when there is a near mass degeneracy with the LSP since there is a strong Boltzmann 
suppression factor Bj in Eq. |^ . 


Ki{{mi + mj)/Tf) 
Ki{2m^o/Tf) 



(mj+TTij — 2m, ..Q ) 
^1 

m ,n 


( 2 . 11 ) 


where rrii^mj are the masses of the incoming particles. To speed up the program a given 
subprocess is removed from the sum (|2.6|) if the total mass of the incoming particles is such 
that Bf is below some limit, B^ dehned by the user. In most cases the cross-section for 
a subprocess with coannihilation is of the same order as the one for the main subprocess 
with only the LSP, acoan ~ Sive a precision of 1% it would be sufficient to 

use = 10“^. However, in our calculation we use a more restrictive value, = 10“®, 
to be safe and to allow for cases where acoan S> occur for example for 

coannihilation processes with squarks which depend on as, for processes with poles or 
in some regions of parameter space where cr^o^o is suppressed. This criteria, in the case 
Xf = 25, corresponds to a restriction on the masses of the incoming particles 


rrii + rrij < 2.55 m^o 


( 2 . 12 ) 


In our program we provide two options to do the integrations, the fast one and the 
accurate one. The fast mode already gives a precision of about 1% which is good enough 













for all practical purposes. The accurate mode should be used only for some checks. 
In the accurate mode the program evaluates all integrals by means of an adaptative 
Simpson program. It automatically detects all singularities of the integrands and checks 
the precision. In the case of the fast mode the accuracy is not checked. We integrate the 
squared matrix elements over cos 6*, the scattering angle, by means of a 5 point Gauss 
formula. For integration over s, Eq. we use a restricted set of points which depends 
whether we are in the vicinity of a s-channel IIiggs(Z,W) resonance or not. We increase 
the number of poins if the Boltzmann factor corresponding to rUpoie is larger than 0.01 

2.3 Higgs width 

When the LSP is near a Higgs resonance, they annihilate very efficiently. The value of 
the neutralino annihilation cross-section depends on the total width which can be very 
large at Higgs masses of ITeV or so. This is particularly so at large tan/3 due to the 
enhancement in the bb channel. However the width of h{H, A) —>■ bb receives important 
QCD corrections. Typically for the heavy Higgses(> ITeV) the partial width into qq can 
vary easily by a factor of 2 from the tree-level prediction, due mostly to the running of the 
quark mass at high scale. To take these corrections into account we have redehned the 
vertices hqq, Hqq and Aqq using an effective mass that reproduces the radiatively corrected 
width. We used HDECAY|^ to produce a table of mass-dependent QCD-corrected Higgs 
partial widths.^ From this, the effective quark masses are extracted and simple 

interpolation can reproduce the one-loop corrected width for any value of rriH- In the 
region of physical interest, the precision on the width is at the per-mil level safe for 
the neutral Higgses partial widths near the ft threshold. However, this region does not 
contribute signihcantly to the neutralino cross-section. For the charged Higgs, one can 
extract from HDECAY both an effective rrit as well as an effective rub using high and low 
tanjS values. This way we could reproduce at better than 1% level the partial width 
—>■ tb. The effective quark mass mb{Q) evaluated at Q = mi + m 2 , the sum of the 
initial masses, is used by default in the {h, H, A) — bb vertices. This will lead to the correct 
result when the contribution of the Higgs resonance is very important {Mh ^ 2m^o). 

2.4 Convention for MSSM parameters 

The input parameters are the ones of the soft SUSY Lagrangian dehned at the weak 
scale, using the same notation as in CompHEP/SUSY modelsj^. In the model used, the 
masses of fermions of the hrst generation are set to zero. The masses of the quarks of 

^We have neglected the SUSY decay modes of the Higgses. These are small in the cases where the 
contribution of the Higgs channel to the relic density is very important, that is when mn ~ 



the second generation are also set to zero, so that there is no mixing of the squarks 
of the first two generations. However we have kept the mass for the muon as well as 
the trilinear coupling as this is relevant for the calculation of the muon anomalous 
magnetic moment. Although hrst generation sleptons are pure left/right states, they are 
properly ordered according to their masses so that the correct coannihilating particle 
corresponding to the lightest slepton is taken into account. In order to avoid unnecessary 
confusion in the sign convention for the parameters fi and A, we give explicitly the mass 
matrices for SUSY particles. 

The physical masses for charginos and neutralinos are obtained after diagonalizing the 
mass matrices, for charginos: 


■Mr. — 


M 2 y/2Mw sin (5 

\/2Mw cos (3 p 


(2.13) 


and for neutralinos 
/ 

Af At = 


Ml 0 —Mz cos (3sw Mz sin (3sw ^ 

0 M 2 Mz cos j3cw —Mz sin j3cw 

Mz cos Psw Mz cos jScw 0 —/i 

V MzAnjSsw —Mz sin (3cw —/i 0 / 


(2.14) 


For sfermions, the mixing becomes signihcant only for the third generation, for t- 
squark, the mass matrix reads 


2_/ + MzCos2p{^- 


Mi = 


mt{At — /i cot f3) 


mt{At - fi cot(3) 
m? +m‘f + |M| cos 2f3s‘^ 


(2.15) 


while for b-squarks 


Mi = 


and staus. 


cos 2/9(-i + 2s^) 
mh{Ah — /itan/?) 


mb{Ab — fitaii(3) 

+ ml — |M| cos 2/?s^ 


J ^2 _ f + M| COS 2/9(-| + s^) mr(Ar - tan/3) 

^ y mriAr —/u tan/3) — M| cos2/ds^ 


(2.16) 


(2.17) 


3 Contents of the micrOMEGAs Package 

micrOMEGAs is a C program that also calls some external FORTRAN functions. micrOMEGAs 
relies on CompHEPf^ for the definition of the parameters and the evaluation of all cross- 
sections. In the package three sample main programs are provided. The hrst uses an 




input file to read the weak scale soft SUSY parameters, the last two give the user the 
option of working in the context of a SUGRA-type model. In the latter case, a sample 
hie containing a call to ISASUGRAQ to run the renormalization group equations from the 
GUT scale to the weak scale is provided. This assumes that the user has already in¬ 
stalled ISASUGRA independently [^. The user has the choice to substitute his/her own 
favourite RGE code. micrOMEGAs is a self-contained program with all subroutines in¬ 
cluded safe for the ISASUGRA routines. The package includes routines for the input of the 
MSSM/SUGRA parameters, calls to CompHEP to calculate physical parameters, calcula¬ 
tion of as well as calls to other subroutines that check for other constraints on the 
model such as 6 —> S 7,((7 — 2 )^ as well as direct limits on the masses of the sparticles. 


Table 1: Weak scale input parameters in the notation of CompHEP 
mu fi Higgs mass parameter 
MGl Ml Gaugino mass 
MG2 M 2 Gaugino mass 
MGS M 3 Gaugino mass 

Mil Left-handed slepton mass, 1st generation, 

M12 Left-handed slepton mass, 2nd generation, 

MIS Left-handed slepton mass, Srd generation, 

Mrl Right-handed selectron mass, 

Mr2 Right-handed smuon mass, 

MrS Right-handed stau mass, 

Mql Left-handed squark mass, 1st generation, 

Mq2 Left-handed squark mass, 1st generation, 

MqS Left-handed squark mass, 1st generation, 

Mul Right-handed u-squark mass, 

Mu2 Right-handed c-squark mass, rric,^ 

MuS Right-handed t-squark mass, 

Mdl Right-handed d-squark mass, 1st generationmj^ 

Md2 Right-handed s-squark mass, 

Md3 Right-handed b-squark mass, 

Atop Trilinear coupling for top quark. At 
Ab Trilinear coupling for b-quark, Ab 
Atau Trilinear coupling for r-lepton, Ar 
Am Trilinear coupling for /i, 

MH3 Mass of Pseudoscalar Higgs {Ma) 

tb Ratio of the vacuum expectation values of the scalar doublets 


Some test input hies with the MSSM weak scale parameters are stored in the data/ 
and datap/ directories together with the reference output hie for comparison purposes. 


"‘The link to Isajet is done here for Version 7.58. 






3.1 Input parameters 


Weak scale parameters 


When running micrOMEGAs, the input hie must specify the set of weak scale SUSY 
parameters as in Table |^. The sfermion masses correspond to the soft masses without 
the D-term. The convention for the third generation of sfermions is specihed explicitly 
in Eq. |2.15| . As in the CompHEP/SUSY model[^|, only the mass matrices for the third 
generation will then be diagonalized to obtain the physical masses. 


Table 2: Standard Model Parameters 


EE 

0.31223 

Electromagnetic coupling constant 

GG 

1.238 

Strong coupling constant at 

SW 

0.473 

sin of the Weinberg angle, sinOw 

sl2 

0.221 

sin of Gabibbo angle 

MZ 

91.1884 

Z mass 

Mm 

0.1057 

mass of muon 

Mtau 

1.777 

mass of tan-lepton 

Me 

1.42 

c-quark mass 

Mtop 

175 

top quark mass 

Mb 

4.62 

b-quark mass 

wtop 

1.7524 

width of top quark 

wZ 

2.4944 

Z-boson width 

wW 

2.08895 

W-boson width 


In addition to these parameters one can choose to redehne parameters, such as the 
standard model parameters (Table H), that are hxed by default. One can also redehne 
SUSY parameters that do not appear in Table 1, such as the widths of Higgses or SUSY 
particles. For example, in processes with t-channel poles it is sometimes necessary to 
specify a width for some particles such as gauginos. By default these widths have been 
set to IGeV. Even though the LSP is assumed to be stable, to avoid any spurious pole it 
is necessary to introduce also a small width. Numerically, the results do not depend on 
the exact value chosen for these widths. A complete list of widths for all supersymmetric 
particles can be found in Table 

SUGRA GUT-scale Parameters 

The value of the weak scale soft supersymmetric Lagrangian can be extracted from 
ISASUGRA. To use this option it is necessary to specify the path for the Isajet libraries. 
The only difference between this option and the weak scale parameters we have just 
mentionned is that the readInitFile has been replaced and an additional call to the 
sugraVar function is done (see the description of these functions in Section 3.2). 










Table 3: Names, masses and widths of supersymmetric particles in CompHEP/SUSY notation 


Name 

symbol 

mass 

width 

Name 

symbol 

mass 

width 

chargino 1 

'1+ 

MCI 

wCl 

mu-sneutrino 

'nm 

MSnmu 

wSnmu 

chargino 2 

'2+ 

MC2 

wC2 

tau-sneutrino 

'nl 

MSntau 

wSntau 

neutralino 1 

~ol 

MNEl 

wNEl 

u-squark L 

'uL 

MSuL 

wSul 

neutralino 2 

~o2 

MNE2 

wNE2 

u-squark R 

'uR 

MSuR 

wSu2 

neutralino 3 

~o3 

MNE3 

wNE3 

c-squark 1 

'cl 

MScl 

wScl 

neutralino 4 

~o4 

MNE4 

wNE4 

c-squark 2 

'c2 

MSc2 

wSc2 

gluino 

~g 

MSG 

wSG 

t-squark 1 

'tl 

MStopl 

wStopl 

Solectron 1 

~el 

MSel 

wSel 

t-squark 2 

't2 

MStop2 

wStop2 

Solectron 2 

'e2 

MSe2 

wSe2 

d-squark L 

'dR 

MSdL 

wSdl 

smuon 1 

~ml 

MSmul 

wSmul 

d-squark R 

'dL 

MSdR 

wSd2 

smuon 2 

~m2 

MSmu2 

wSmu2 

s-squark L 

'sL 

MSsL 

wSsl 

stau 1 

'll 

MStaul 

wStaul 

s-squark R 

'sR 

MSsR 

wSs2 

stau 2 

'12 

MStau2 

wStau2 

b-squark 1 

'bl 

MSbotl 

wSbotl 

e-sneutrino 

'ne 

MSne 

wSne 

b-squark 2 

'b2 

MSbot2 

wSbot2 


Table 4: Parameters for Higgs particles in CompHEP/SUSY notation 


Name 

symbol 

mass 

width 

Name 

symbol 

mass 

width 

Light Higgs 

h 

Mh 

wh 

CP-odd Higgs 

H3 

MH3 

wH3 

Heavy higgs 

H 

MHH 

wHh 

Charged Higgs 

H+ 

MHc 

wHc 


The set of input parameters for the SUGRA model are the usual 

mO common scalar mass at GUT scale 

mhf common gaugino mass at GUT scale 

aO trilinear soft breaking parameter at GUT scale 

tb tan(beta) or the ratio of vacuum expectation values 

sgn +/~1, sign of Higgsino mass term 

mtop pole mass of the top quark 

model = 1 SUGRA model 

= 2 SUGRA with true gauge coupling unification at the GUT scale 

3.2 Functions of micrOMEGAs 

All functions provided with the package are briefly described below. The user can hnd 
examples of how to call these routines in the omg. c or sugomg. c hie located in the root 
directory. The types of functions called are specihed in Table 

• darkOmega(Xf,fast,Beps) 

This is the basic function of the package which returns the relic density VLh? (Eq. |2.3|) . 
Xf is the initial value of the Xf parameter needed for solving iteratively the freeze-out 






















Table 5: Main functions called 


double darkOmega(double fast, double Beps) 

int readlnitFile(char * filename) 

int assignVal(char * name, double val) 

int sugraVar(float mO,float mhf,float aO,float tb,float sgn, float mtop, 
int model) 

int calcDep(int LC, int Widths) 
int findVaKchar * name, double * val) 
void printMasses(FILE * f, int sort) 
char * Isp(void) 

void printMainChannels(FlLE *f, double cut) 

void printMSSMvar(FILE *f) 

double delrho(void) 

double gmuon(void) 

int MassLimits(void) 

double bsganew(void) 


temperature, Eq. ( |2.8|) . We recommend to use Xf = 25 for the first call, this is the value set 
by default in the sample main hie. The value of Xf returned by darkOmega corresponds to 
the solution of the freeze-out equation (|2.8| ). If the user calls darkOmega in a scan where 
the model parameters are modihed slowly, we recommend to use the preceding Xf value 
as input for the next iteration. Otherwise, it is better to assign to Xf some hxed value 
before each call. The parameter Beps dehnes the criteria for including a given channel 


m 


the sum for the calculation of the thermally averaged cross-section, Eq. (|2.11|) , 10 


-6 


IS 


the recommended value. 


If fast 7 ^ 0, the program uses the fast mode, otherwise the accurate mode is used, 
see Section 2.2. If some problem is encountered, darkOmega returns -1. 

• readInitFile(filename) 

Reads the input hie, filename. The input parameters of the MSSM that need to be 
specihed are given in Table |^. The function returns: 

0 - then the input has been read correctly; 

-1 - if the hie does not exist or can not be opened for reading; 
n - when the line number n has been written in the wrong format. 

The correct format of a line is 

<name of variableXspaceXnumerical value> 
for example 

MGS 1500. 

This function writes a WARNING if one of the variables of Table 1 is not dehned in 
filename 







•assignVal(name,val) 


Changes the value of one of the standard model or MSSM parameters, name is either 
one of the parameter names given in Tables or a width of a particle whose name is 
specified in Tables val is the new value. The function returns 0 when it successfully 
recognizes the parameter name and 1 otherwise. The only reason for an error is a wrong 
parameter name. 


•calcDep(LC,Widths) 


Calculates all parameters of the MSSM model including physical masses and mixings as 
well as couplings. The calculation of the Higgs masses is based on FeynHiggsf astpO[ . 


calcDep returns a non zero error code when it is unable to calculate some physical 
parameter. For example, this can occur when the input parameters are such that one ends 
up with a negative squared mass for a sfermion. The input parameters of calcDep should 
be set to LC=1 and Width=l. We provide this option only to give the possibility to make 
comparisons with other programs. LC=0 switches off the QCD loop correction for the 
Higgs bosons decays to b, c,and t quarks. Width=0 switches off the automatic calculation 
of widths of Higgs bosons. If Width=0 the user has to define the Higgs width values by 
means of either readInitFile or assignVal. The names of the corresponding variables 
are listed in Table Note that the choice for the Width parameter affects only the Higgs 
width in the propagator. The value of the Hqq vertex is calculated independently either 
at the tree-level(LC=0) or one-loop level(LC=l) and could be different from the one that 
reproduces the Higgs total width. 


We have to emphasize that the user must call calcDep after having defined the free 
MSSM parameters and before the calculation of the relic density. Otherwise the calculated 
relic density is not valid. 


•findVal(name,feval) 


Outputs the value of a parameter of the MSSM, either one of the input parameters or one 
of the derived parameters, name should be the name of a variable (free or constrained one). 
findVal returns zero if the name indeed corresponds to some variable and 1 otherwise. 
For a name that has been correctly specified, f indVal assigns to val the numerical value 
of the parameter. This function is particularly useful if one wants to impose constraints on 
the masses of the superparticles. The names of the masses of superparticles and Higgses 
together with the names of sparticles in CompHEP notation can be found in Tables |. 

•printMasses(filename, sort) 

^Note that some of the mass parameters can be negative, the masses of particles are automatically 
redefined such that they correspond to the absolute values for these parameters. 





This routine prints the masses of the supersymmetric particles mentionned in Table ^ 
as well as all Higgs masses and widths into the specified file filename. To see the masses 
on the screen replace filename by stdout. If sortT^ 0 the masses are sorted in order of 
increasing mass. 

•printMSSMvar(filename) 

Prints the MSSM parameters presented in the Table 1. To see the parameters on the 
screen replace filename by stdout. 

•printMainChannels(filename,cut) 

Allows the user to see which channels dominate in the calculation of VLh?. This function 
writes into the file filename, the channels whose relative contribution to l/Hh^ exceeds 
the value chosen for cut. The value of the contribution in percent is also written. If 
cut=0, all channels will be displayed. 


•sugraVar(mO,mhf,aO,tb,sgn,mtop,model) 

This function solves the renormalization group equations and assigns a value to all the 
MSSM weak scale parameters in Table 1, except for the trilinear coupling for the muon, 
Am (Afj). We have used the approximation, = Aq — .7 * M 1/2 to dehne this parameter. 
This function calls the Isajet package, which must be installed by the user. sugraVar 
returns zero for a succesful assignment of weak scale parameters or the Isajet error code 
for a wrong set of parameters. The error code corresponds to the NOGOOD parameter of 
COMMON /SUGPAS/, (see the Isajet manual |]T9|). To substitute another RGE code, 
the user must modify the files sources/sugrac.c. 


•IspO 

This function returns the name of the LSP. The relic density can be calculated with 
any particle being the LSP. If the user wants to impose a specific LSP, the nature of the 
LSP must be checked after calling calcDep. In the sample omg.c hie that is provided in 
the package, this is done for the conventional choice of a neutralino LSP. 


•bsganewO 


Returns the value of the branching ratio for b —>■ sj. For 6 —> 57 we have improved 
on the results of by including some very recent new contributions beyond the leading 
order that are especially important for high tan/3. We include these terms by following 
This function should be called after all parameters of the MSSM have been dehned 


in calcDep. 


•gmuonO 


Returns the value of the supersymmetric contribution to the anomalous magnetic 




moment of the muonpQt ph| . The result depends only on the parameters of the gaugino 
mass matrices as well as on the smuon parameters, in particular the trilinear coupling. 
This parameter has been added in the function sugraVar described above as this is not 
calculated by ISASUGRA/Isajet. 

•delrhoO 


Computes the stop/sbottom contribution to the Ap parameter, including the two- 
loop gluon exchange and pure scalar diagrams]^. Returns the value for Ap and a 
WARNING if the value exceeds the limit Ap > 1.3 x 10“^. This routine is included 
in FeynHiggsFast pPj] . 

•MassLimitsO 


This function returns a value greater than 1 and prints a WARNING when the choice 
of parameters conflicts with a direct accelerator limits on sparticle masses. Only the 
direct limits from LEP in the general MSSM are implemented. In constrained models 
more severe limits can be set, this should be added by the user. The constraint on the 
light Higgs mass is not implemented and must be added by the user. 


4 Compilation and installation instructions. 

The package can be obtained at 

http://wwwiapp.in2p3.fr/lapth/micromegas 

After unpacking the hies, the root directory of the package, 
micromegas.l.0 

will be created. This directory contains a makehle, three sample main programs omg. c, 
sugomg.c and scycle.c and two directories for data hies. The compiler options and 
the location of the Isajet libraries, when needed, should be dehned in the makehle. 
The compiler options for Linux, Alpha and Silicon Graphics workstations are dehned 
automatically. At the moment the package cannot be compiled on HP Workstations. 

To compile, type 
make 

This creates the omg executable hie for the omg.c main program which does not require 
the Isajet library. If omg is launched without an argument it calculates the relic density 
for the default parameters which corresponds to the data/datapl hie MSSM. The output 
is presented in Section 5. 

If omg is launched with some arguments, it treats all of them as names of input hies 







and calculates the relic densities for each of them. Say, 
omg datap/* 

produces the results presented in the Table 

As the compilation of the code for more than 2800 processes would take a long time 
while in most cases only a small fraction of processes are needed, we have generated 
the processes and linked them dynamically only when they are needed. The hrst time 
new processes need to be linked compilation can take up to a few minutes, the actual 
computation time remains insignihcant. 

The sugomg.c code contains an example of application of the sugraVar routine. It 
can be compiled by 

make main=sugomg 

The corresponding executable hie is sugomg. It is necessary to hrst install the Isajet pack¬ 
age. The environment variable SUGRA is used to dehne the path to the isaj et. a library. 
This variable should be dehned properly by modifying the makefile. 

A sample hie scycle.c for making a scan over two parameters (mo,mi/ 2 ) in the 
SUGRA model is also provided. The range of values for mo,mi /2 need to be specihed, 
the remaining parameters have to be given as inputs. It can be compiled by 
make main=scycle 

The corresponding executable hie is scycle. 

To install the package one needs initially about 10MB of disk space. As the program 
generates libraries for annihilation processes only at the time they are required, the total 
disk space necessary can increase signihcantly after running the program for diherent 
models. These libraries are stored in the sources/2-2/ subdirectory and can be removed 
by the user. However, this is done at the expense of slowing down the compilation 
whenever the program will need again these libraries. 

5 Sample Output 

Running micrOMEGAs with the default value of the standard parameters and with the 
input hie datap/datapl, which contains the parameters of the MSSM as specihed in 
Table model A, will produce the following output: 

Initial file "data/datapl" 


Higgs masses and widths 
h : Mh =125.1 (wh 


=3.2E-03) 


H 

: MHH 

= 1000.8 

(wHh 

=2.3E+00) 

H3 

: MH3 

= 1000.0 

(wH3 

=2.4E+00) 

H+ 

: MHc 

= 1003.1 

(wHc 

=2.3E+00) 


Masses of SuperParticles: 


~ol 

: MNEl = 

22.4 1 

~ne 

: MSne 

80.0 

1 1 

~nm 

: MSnmu = 

80.0 

~11 

: MStaul= 

92.7 1 

~nl 

: MSntau= 

93.1 

1 1 

~ml 

: MSmul = 

108.0 

~el 

: MSel = 

108.2 1 

~e2 

: MSe2 = 

111.2 

1 1 

~m2 

: MSmu2 = 

111.3 

~12 

: MStau2= 

124.7 1 

~1+ 

: MCI 

198.6 

1 1 

~o2 

: MNE2 = 

200.1 

~o3 

: MNE3 = 

270.3 1 

~g 

: MSG 

300.0 

1 1 

~2+ 

: MC2 

332.4 

~o4 

: MNE4 = 

332.6 1 

~tl 

: MStopl= 

813.0 

1 1 

~bl 

: MSbotl= 

998.4 

~cl 

: MScl = 

998.7 1 

~uL 

: MSuL = 

998.7 

1 1 

~uR 

: MSuR = 

999.4 

~c2 

: MSc2 = 

999.4 1 

~dR 

: MSdR = 

1000.3 

1 1 

~sR 

: MSsR = 

1000.3 

~dL 

: MSdL = 

1001.6 1 

~sL 

: MSsL = 

1001.6 

1 1 

~b2 

: MSbot2= 

1003.5 

~t2 

: MStop2= 

1175.8 1 









Channels which contribute more than 1 °/o. 

257o ~ol ~ol -> e E 
257 ~ol ~ol -> m M 
397 ~ol ~ol -> 1 L 
47 ~ol ~ol -> ne Ne 
47 ~ol ~ol -> nm Nm 
27 ~ol ~ol -> nl N1 
0mega=2.48E-01 
g2=2.55E-09 
bsgamma=4.4E-04 

Running micrOMEGAs in the SUGRA case, with ISASUGRA version 7.58 and the input 
parameters 

sugomg 100 400 0 10 +1 175 1 

will write the values of the input parameters, the value of the weak scale parameters 
as calculated by ISASUGRA, the corresponding list of weak scale parameters in CompHEP 
notation followed by the usual output of micrOMEGAs as described above: 

MSSM parameters: 
mu 513.940063 

MGl 162.922958 











MG2 311.489258 

MGS 914.497253 

Mil 286.277618 

M12 286.277618 

M13 285.412537 

Mrl 178.365540 

Mr2 178.365540 

Mr3 175.411743 

Atau -244.407349 

Am -280.000000 

Mql 819.445862 

Mq2 819.445862 

Mq3 760.147278 

Mul 793.780884 

Mu2 793.780884 

Mu3 659.573975 

Mdl 790.153137 

Md2 790.153137 

Md3 785.846619 

Atop -723.817383 

Ab -1038.150269 

MH3 576.585266 


Higgs masses and widths 


h 

: Mh 

= 115.2 

(wh 

=3.2E-03) 

H 

: MHH 

= 576.9 

(wHh 

=1.2E+00) 

H3 

: MH3 

= 576.6 

(wH3 

=1.3E+00) 

H+ 

: MHc 

= 582.2 

(wHc 

=1.3E+00) 


Masses of SuperParticles: 


~ol : MNEl = 160.8 

~el : MSel = 183.4 

~nm : MSnmu = 279.1 

~12 : MStau2= 292.3 

~o3 : MNE3 = 518.2 

~tl : MStopl= 623.2 

~sR : MSsR = 790.5 

~b2 


~11 

: MStaul= 

175.7 

~nl 

: MSntau= 

278.2 

~e2 

: MSe2 = 

290.2 

~1+ 

: MCI 

297.0 

~o4 

: MNE4 = 

534.5 

~bl 

: MSbotl= 

752.2 

~uR 

: MSuR = 

793.0 

~uL 

: MSuL = 

817.7 


~ml 

: MSmul = 

183.4 

~ne 

: MSne 

279.1 

~m2 

: MSmu2 = 

290.2 

~o2 

: MNE2 = 

297.3 

~2+ 

: MC2 

534.6 

~dR 

: MSdR = 

790.5 

~cl 

: MScl = 

793.0 

~c2 

: MSc2 = 

817.7 


MSbot2= 796.1 











~t2 : MStop2= 818.3 II ~dL : MSdL = 821.6 II ~sL : MSsL = 821.6 

~g : MSG = 914.5 II 


Channels 

which ( 

contribute more than 1 

17. 

~ol 

~ol 

-> 

b 

B 

177. 

~ol 

~ol 

-> 

e 

E 

177. 

~ol 

~ol 

-> 

m 

M 

187. 

~ol 

~ol 

-> 

1 

L 

57. 

~ol 

~11 

-> 

Z 

1 

207. 

~ol 

~11 

-> 

A 

1 

17. 

~ol 

~ml 

-> 

Z 

m 

57. 

~ol 

~ml 

-> 

A 

m 

17. 

~ol 

~el 

-> 

Z 

e 

47. 

~ol 

~el 

-> 

A 

e 

37. 

~11 

~11 

-> 

1 

1 

0mega=2.07E-01 



g2=l. 

,46E- 

-09 





bsgaiiima=3.48E-04 


6 Results and Comparisons 


The micrOMEGAs code was extensively tested against another pnblic package for calcn- 
lating the relic density, DarkSUSY. As discnssed previously, the two codes differ somewhat 
in the numerical method used for solving the density equations. micrOMEGAs includes 
more subprocesses(e.g. all coannihilations with sfermions), loop-corrected Higgs widths, 
and complete tree-level matrix elements for all processes. DarkSUSY includes on the other 
hand some loop induced processes such as yy —> gg, 77 which are generally small. When¬ 
ever the coannihilation channels with sfermions and the Higgs pole are not important 
we expect good agreement with DarkSUSY. We have first compared micrOMEGAs with 
a version of DarkSUSY where we have replaced the matrix elements by CompHEP matrix 
elements. As expected, we found excellent numerical agreement between the two pro¬ 
grams (at the 2% level) for all points tested. We have then made comparisons with the 
original version of DarkSUSY. The results of these comparisons are displayed for a few test 
points in Table ^ We found in general good agreement between the two codes. However 
we have observed some discrepancies that could reach up to (30%) in particular in the 
process XiXi t&(model B in Table ^ H As displayed in the line flfree^ when removing 


“The CompHEP result for this matrix element agrees with the result of GraceSUSY|31 . 




non-gaugino coannihilation channels and reverting to the tree-level treatment of the Higgs 
width we recover results similar to DarkSUSY. The impact of these extra channels, model 
C for sleptons and model D for squarks can be as large as an order of magnitude and 
depends critically on the mass difference with the lightest neutralino. 


Table 6; Sample results and comparison with DarkSusy 


Model 

A 

B 

G 

D 

E 

F 

tb 

5. 

10 . 

10 . 

10 . 

45. 

50. 

mu 

264.5 

400. 

518.6 

-1200. 

500. 

1800. 

MGl 

25.9 

500. 

166.1 

300. 

180. 

850. 

MG2 

258.9 

1000 . 

317.9 

600. 

350. 

1600. 

MG3 

800. 

3000. 

931.8 

1800. 

1000 . 

4000. 

Mil 

101 . 

1000 . 

289.0 

1000 . 

500. 

2000 . 

M13 

112 . 

1000 . 

288.1 

1000 . 

500. 

2000 . 

Mrl 

100 . 

1000 . 

177.0 

1000 . 

300. 

1600. 

Mr3 

88 . 

1000 . 

174.1 

1000 . 

250. 

1600. 

Mql 

1000 . 

1000 . 

834.1 

1000 . 

1000 . 

3500. 

Mq3 

1000 . 

1000 . 

773.9 

1000 . 

1000 . 

3500. 

Mul 

1000 . 

1000 . 

803.1 

500. 

1000 . 

3500. 

Mu3 

1000 . 

1000 . 

671.8 

500. 

1000 . 

3500. 

Mdl 

1000 . 

1000 . 

799.3 

500. 

1000 . 

3500. 

Md3 

1000 . 

1000 . 

799.7 

500. 

1000 . 

3500. 

Atop 

2400. 

0 . 

-738.1 

-1800. 

-1000. 

-3000. 

Ab 

2400. 

0 . 

-1058. 

-1800. 

-1000. 

-3000. 

At an 

0 . 

0 . 

-249.2 

0 . 

-100. 

-500. 

Mh3 

1000 . 

1000 . 

581.9 

1000 . 

500. 

1700. 

m^o 

22.4 

384.3 

164. 

299.9 

178.2 

849.3 


200.3 

402.7 

303.5 

597.9 

333.9 

1585.4 


198.6 

395.5 

303.3 

597.9 

333.8 

1585.3 

rriei 

108.2 

1000.9 

182.1 

1000.9 

303.0 

1600.6 

rrif. 

92.7 

997.5 

174.4 

990.3 

236.9 

1595.0 

mu 

813.0 

1007.7 

635.6 

326.8 

930.7 

3439.9 

rrih 

125.1 

111.2 

115.3 

118.4 

118.1 

121.7 


25.5 

1.88 

13.9 

-1.66 

30.4 

2.07 

b S7(xl0“^) 

4.4 

3.8 

3.5 

4.9 

2.8 

3.5 

micrOMEGAs 

.25 

.024 

.14 

.15 

.22 

.26 

QX 
^ ^tree 

.25 

.024 

.32 

.74 

.13 

.56 

DarkSusy 

.25 

.018 

.32 

.74 

.13 

.55 


The effect of the Higgs width is particularly important at large tan/? with the en¬ 
hanced contribution of the &-quark. However the one-loop correction amounts to a re¬ 
duced effective b-quark mass and a much smaller width especially at large values of mn- 
If it was not for the strong Boltzmann suppression factor singling out the contribution 
at y/s ~ 2m^ there will be little difference after integrating over the peak for the one- 


loop or tree-level result. However the effect observed can be as much as a factor 2. For 
~ Ma/2 the narrower resonance suffers less from the Boltzmann reduction factor 
leading to < cr^-’-°°P > / < atree >> 1 and < VLtree (model F). Further away from 

the pole however one catches the contribution from the wider resonance without excessive 
damping from the Boltzmann factor, as expected Vti_ioop > ^tree (model E). 


Our numerical results were also compared with Ref. [^. Qualitative agreement 


is found in the case of SUGRA model although we use a different RGE code. Precise 
comparisons necessitates a careful tuning of parameters to make sure we have the same 
parameters at the weak scale. A random scan over mo — mi /2 for /r > 0 shows the typical 
shape of the allowed region (.1 < < .3) in SUGRA model for moderate values of 

tan (3. 


Figure 1: VLh? in the mo — ■mi /2 plane in a SUGRA model with tan/? = 10, /r > 0 and 
mtop = 175GeV. The dark (blue) region corresponds to .1 < VLh? < .3 and the light 
(pink) region to VLh? < .1. In the white region at large mo — Tni/ 2 : the relic density is 
above the present limit G > .3 and in the region at small tuq, the f is the LSP. The 
yellow (light grey) region correspond to the region excluded by 6 —sy. The LEP limit 
on mh = 113GeV is also displayed. 



Mi/2(GeV) 









7 Conclusion 


The package micrOMEGAs allows to calculate the relic density of the LSP in the most 
general MSSM with Rp conservation. The package is self-contained safe for the ISASUGRA 
package that is required when using the SUGRA option. All possible channels for coan¬ 
nihilations are included and all matrix elements are calculated exactly at tree level with 
the help of CompHEP. Loop corrections for the masses of Higgs particles (two-loop) and 
the width of the Higgs (QCD one-loop) are implemented. Good agreement with existing 
calculations is found when identical set of channels are included. Future versions will 
include loop corrections to neutralino masses. Even though these corrections are only a 
few GeV’s they can alter signihcantly the calculation of the relic density when there is a 
near mass degeneracy with the next to lightest supersymmetric particle that contributes 
to a coannihilation channel . Although the loop processes are in general small, we will 
include yy —>• 77 , 7 Z, gg in an update of micrOMEGAs. 
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